#### setting environment ####
require(runjags)
require(coda)
require(lmtest)
require(sandwich)
require(multiwayvcov)
require(dichromat)
color <- colorschemes$Categorical.12  # color for figures
week.id.function <- function(date, end.date) {  # function to identify the week-ID of a date
  for (i in 1:length(end.date)) {
    if (date <= end.date[i]) {
      return(i)
      break
    }
  }
}
posterior.summary <- function(x) {  # function to compute a posterior mean and a 95% HPD interval
  posterior.mean <- mean(x)
  HPD.interval <- HPDinterval(as.mcmc(x))
  c(posterior.mean, HPD.interval)
}
coef.table <- function(sample.matrix) {  # function to create a coefficient table
  summary.table <- matrix(NA, 2 * ncol(sample.matrix), 1)
  rownames(summary.table) <- rep("", 2 * ncol(sample.matrix))
  for (i in 1:ncol(sample.matrix)) {
    summary.table[2 * i - 1, 1] <- sprintf("%5.3f", mean(sample.matrix[, i]))
    summary.table[2 * i, 1] <- paste("(", sprintf("%5.3f", sd(sample.matrix[, i])), ")", sep = "")
    rownames(summary.table)[2 * i - 1] <- colnames(sample.matrix)[i]
  }
  summary.table
}

#### preparing the data and variables ####
## loading the polling data
setwd("weekly_approval_rating")
# confirm that the current working directory contains only 11 polling data files
filename <- dir()
week.id.data <- read.csv(filename[1])[, c(1:3)]
week.id.data$begin <- as.Date(week.id.data$begin)
week.id.data$end <- as.Date(week.id.data$end)

weekly.data <- c()
for(i in 1:(length(filename))){
  weekly.data <- rbind(weekly.data, read.csv(filename[i]))
}
setwd("../")

## cleaning up the dataset
weekly.data <- na.omit(weekly.data)
weekly.data <- weekly.data[order(weekly.data$pmid), ]  # data is sorted by PMs
rownames(weekly.data) <- NULL

## Table 1 (polling result statistics by polling firms)
table.polling.firms <- matrix(NA, 12, 4)
for (i in 1:11) {
  subdata <- subset(weekly.data, house == i)
  table.polling.firms[i, 1] <- nrow(subdata)  # total number of polls
  table.polling.firms[i, 2] <- sum(subdata$reshuffle)  # number of polls immediately following a reshuffle
  table.polling.firms[i, 3] <- sum(subdata$cue == 1)  # number of polls providing a type I reshuffle cue
  table.polling.firms[i, 4] <- sum(subdata$cue == 2)  # number of polls providing a type II reshuffle cue
}
table.polling.firms[12, ] <- colSums(table.polling.firms[1:11, ])
table.polling.firms

## recording where each PM starts and ends in the dataset
start.week <- end.week <- rep(NA, 8)
counter <- 0
for (p in 1:8) {
  start.week[p] <- counter + 1
  end.week[p] <- counter <- counter + sum(weekly.data$pmid == (p + 55))
}

## the upper and lower limits of the prior distribution of each PM's approval rating in the first week
# the result of Jiji Press's opinion poll conducted for the first time after each PM was inaugurated plus/minus 0.2
lo <- up <- rep(NA, 8)
for(p in 1:8){
  weekly.pm <- subset(weekly.data, pmid == (p + 55) & house == 1)
  lo[p] <- max(0, (weekly.pm$result[weekly.pm$weekid == min(weekly.pm$weekid)] - 20) / 100)
  up[p] <- min(1, (weekly.pm$result[weekly.pm$weekid == min(weekly.pm$weekid)] + 20) / 100)
}

## creating variables for the main analysis
# ID of the week when each PM was inaugurated
inauguration.week <- c(1, 284, 336, 388, 439, 477, 541, 610, nrow(week.id.data) + 1)
# recoding the ID of weeks so as to the week when each PM was inaugurated takes 1
weekid <- weekly.data$weekid - inauguration.week[(weekly.data$pmid - 55)] + 1
# number of weeks in which each PM was in office
n.weeks <- rep(NA, 8)
for (p in 1:8) {
  n.weeks[p] <- inauguration.week[p + 1] - inauguration.week[p]
}
# date of reshuffles
reshuffle.date <- as.Date(c("2002/9/30", "2003/9/22", "2004/9/27", "2005/10/31", "2007/8/27", 
                            "2008/8/2", "2010/9/17", "2011/1/14", "2012/1/13", "2012/6/4", 
                            "2012/10/1", "2014/9/3", "2015/10/7"))
# week-ID of reshuffles
reshuffle.week.id <- rep(NA, length(reshuffle.date))
for (i in 1:length(reshuffle.date)) {
  reshuffle.week.id[i] <- week.id.function(reshuffle.date[i], week.id.data$end)
}
# indicator variable for in which week each PM reshuffled his cabinet
# for example, "reshuffle[1, 76] <- 1" means Koizumi reshuffled his cabinet in the 76th week
reshuffle <- matrix(0, 8, max(n.weeks))  # the row is for PMs, the column is for weeks
for (p in 1:8) {
  for (i in 1:length(reshuffle.week.id)) {
    if (inauguration.week[p] <= reshuffle.week.id[i] & reshuffle.week.id[i] < inauguration.week[p + 1]) {
      reshuffle[p, reshuffle.week.id[i] - inauguration.week[p] + 1] <- 1
    }
  }
}
# approval rating in each poll
y <- weekly.data$result / 100
# number of respondents in each poll
n.obs <- weekly.data$nobs
# survey mode of each poll
survey.mode <- ifelse(weekly.data$mode == 1, 1,  # face-to-face
                      2)  # telephone
# whether each poll was conducted by RDD or quota sampling
design <- ifelse(weekly.data$mode > 2, 1, 0)
# polling firm of each poll
house <- weekly.data$house
# whether each poll provided a reshuffle cue
cue <- ifelse(weekly.data$cue > 0, 1, 0)

#### estimation of the main analysis ####
## initial values
inits.main <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.1)
  g <- runif(1, -0.03, 0.03)
  d <- runif(10, -0.03, 0.03)
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(1, -0.03, 0.03)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, gamma.star = g, delta.star = d, tau = t, kappa = k, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.main <- list()
for (i in 1:3) {
  initial.values.main[[i]] <- inits.main(1000 + 100 * i, 2000 + 100 * i, i)
}

## MCMC estimation
result.main <- run.jags(model = "JAGS_main.r", 
                        monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                    "alpha5", "alpha6", "alpha7", "alpha8", 
                                    "lambda", "gamma", "delta", "kappa", "omega", "tau"), 
                        data = list(y = y, weekid = weekid, n.weeks = n.weeks, n.obs = n.obs, 
                                    survey.mode = survey.mode, design = design, house = house, 
                                    reshuffle = reshuffle, cue = cue, lo = lo, up = up, 
                                    start.week = start.week, end.week = end.week), 
                        inits = initial.values.main, 
                        n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                        thin = 20, method = "parallel")
# save(result.main, file = "result_main.Rdata")
# load("result_main.Rdata")

# convergence check
which(gelman.diag(result.main$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.main <- rbind(result.main$mcmc[[1]], result.main$mcmc[[2]], result.main$mcmc[[3]])
colnames(sims.matrix.main) <- colnames(result.main$mcmc[[1]])

#### Figure 1 (estimated weekly cabinet approval ratings) ####
## extracting the parameters of approval ratings
approval.samples <- list()
count <- 0
for (p in 1:8) {
  # approval ratings (alpha) are contain in the top of the MCMC sample matrix
  approval.samples[[p]] <- sims.matrix.main[, (count + 1):(count + n.weeks[p])]
  count <- count + n.weeks[p]
}

## computing the posterior means and 95% HPD intervals of weekly approval ratings
ts.plot.data <- list()
count <- 0
for(p in 1:8){
  HPD <- HPDinterval(as.mcmc(approval.samples[[p]]), prob = 0.95)
  ts.plot.data[[p]] <- data.frame(id = (count + 1):(count + n.weeks[p]), 
                                  # posterior mean
                                  rating = colMeans(approval.samples[[p]]), 
                                  # upper and lower bound of the 95% HPD interval
                                  upper = HPD[, 2], lower = HPD[, 1])
  count <- count + n.weeks[p]
}
Jiji.data <- subset(weekly.data, house == 1)  # Jiji's monthly data for comparison

## drawing Figure 1
tiff("Figure_1.tif", width = 150, height = 90, units = "mm", pointsize = 7, res = 600)
par(mar = c(5, 4, 2, 2))
plot(NULL, NULL, type = "n", xlim = c(0, inauguration.week[9] - 1), ylim = c(0, 1), 
     xlab = "", ylab = "Approval Rating", xaxt = "n")
for(p in 1:8){
  # shaded bands representing 95% HPD intervals
  for(i in 1:n.weeks[p]){
    segments(ts.plot.data[[p]]$id[i], ts.plot.data[[p]]$upper[i], 
             ts.plot.data[[p]]$id[i], ts.plot.data[[p]]$lower[i], 
             col = color[11])
  }
  # solod lines representing point estimates
  lines(ts.plot.data[[p]]$id, ts.plot.data[[p]]$rating, col = color[12])
  # Jiji's monthly data for comparison
  sub.Jiji.data <- subset(Jiji.data, pmid == 55 + p)
  lines(sub.Jiji.data$weekid - inauguration.week[1] + 1, sub.Jiji.data$result / 100, 
        lty = 3, col = color[10])
}
# weeks in which PMs reshuffled his cabinet
points(c(76, 127, 180, 237, 332, 380, 491, 508, 560, 581, 598, 698, 755), 
       colMeans(sims.matrix.main[, c(76, 127, 180, 237, 332, 380, 491, 508, 560, 581, 598, 698, 755)]), 
       pch = 4, cex = 0.7)
# horizontal axis
axis(1, at = c(1, 37, 89, 141, 193, 245, 298, 350, 402, 454, 506, 558, 611, 663, 715, 758), 
     labels = c("4/2001", "1/2002", "1/2003", "1/2004", "1/2005", "1/2006", "1/2007", "1/2008", 
                "1/2009", "1/2010", "1/2011", "1/2012", "1/2013", "1/2014", "1/2015", "11/2015"), 
     las = 2)
# PM's labels
text(25, 0.92, labels = "Koizumi")
text(300, 0.73, labels = "Abe (1)")
text(350, 0.62, labels = "Fukuda")
text(400, 0.53, labels = "Aso")
text(450, 0.8, labels = "Hatoyama")
text(490, 0.7, labels = "Kan")
text(550, 0.65, labels = "Noda")
text(620, 0.76, labels = "Abe (2)")
dev.off()

#### Figure 2 (estimated effects of a reshuffle cue, house effects, and mode effects) ####
## extracting parameter samples
lambda.samples <- sims.matrix.main[, 763]
gamma.samples <- sims.matrix.main[, 764:765]
delta.samples <- sims.matrix.main[, 766:776]

## computing the posterior means and 95% HPD intervals
post.lambda.mean <- mean(lambda.samples)
post.lambda.HPD <- HPDinterval(as.mcmc(lambda.samples))

post.gamma.mean <- colMeans(gamma.samples)
post.gamma.HPD <- HPDinterval(as.mcmc(gamma.samples))

post.delta.mean <- colMeans(delta.samples)
post.delta.HPD <- HPDinterval(as.mcmc(delta.samples))

## drawing Figure 2
# data frame for plotting a dotchart
post.plot.matrix <- data.frame(mean = rev(c(post.gamma.mean, post.delta.mean, post.lambda.mean)), 
                               groups = c("Cue Effect", rep("House Effects", 11), rep("Mode Effects", 2)))

tiff("Figure_2.tif", width = 120, height = 150, units = "mm", pointsize = 8, res = 600)
par(mar = c(5, 2, 2, 2))
dotchart(post.plot.matrix$mean, groups = post.plot.matrix$groups, 
         labels = rev(c("Face-to-Face", "Telephone", 
                        "Jiji", "Kyodo", "Yomiuri", "Asahi", "Mainichi", "Sankei-FNN", 
                        "Nikkei", "NHK", "JNN", "ANN", "NNN", "Reshuffle")), 
         pch = 19, xlim = c(-0.07, 0.07), 
         xlab = "Effects (Point)", ylab = "")
points(post.lambda.mean, 18, pch = 19, col = color[12])
points(post.delta.mean, c(15:5), pch = 19, col = color[10])
points(post.gamma.mean, c(2:1), pch = 19, col = color[8])
segments(post.lambda.HPD[1], 18, post.lambda.HPD[2], 18, col = color[12])
segments(post.delta.HPD[, 1], c(15:5), post.delta.HPD[, 2], c(15:5), col = color[10])
segments(post.gamma.HPD[, 1], c(2:1), post.gamma.HPD[, 2], c(2:1), col = color[8])
abline(v = 0, lty = 3)
dev.off()

#### quantity of interests referred in the main text ####
## average change in approval ratings concurrent with reshuffles (kappa)
round(posterior.summary(sims.matrix.main[, "kappa"]), 3)

## effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.main[, "lambda"]), 3)
# how great is the effect compared to the time-series variance of the weekly approval ratings
round(colMeans(sims.matrix.main[, 778:785]), 3)
round(pnorm(mean(lambda.samples), 0, mean(sims.matrix.main[, "omega[6]"])), 2)
round(pnorm(mean(lambda.samples), 0, mean(sims.matrix.main[, "omega[8]"])), 2)

## house effects (delta)
round(post.delta.mean, 3)

## mode effects (gamma)
round(post.gamma.mean, 3)

## Table 4 (estimated results of miscellaneous parameters in the main analysis)
coef.table.4 <- matrix(NA, 19, 3)
for (i in 1:19) {
  coef.table.4[i, ] <- posterior.summary(sims.matrix.main[, 777 + i])
}
rownames(coef.table.4) <- colnames(sims.matrix.main)[778:796]
round(coef.table.4, 3)

#### comparison with the cue of inauguration (online Appnendix B) ####
## variables for the anlaysis of inauguration cues
weekly.data.2 <- subset(weekly.data, house != 9)  # excluding JNN News data
nrow(weekly.data.2)  # number of polls
weekid.2 <- weekly.data.2$weekid - inauguration.week[(weekly.data.2$pmid - 55)] + 1
start.week.2 <- end.week.2 <- rep(NA, 8)
counter <- 0
for (p in 1:8) {
  start.week.2[p] <- counter + 1
  end.week.2[p] <- counter <- counter + sum(weekly.data.2$pmid == (p + 55))
}
y.2 <- weekly.data.2$result / 100
n.obs.2 <- weekly.data.2$nobs
survey.mode.2 <- ifelse(weekly.data$mode == 1, 1, 2)
design.2 <- ifelse(weekly.data.2$mode > 2, 1, 0)
house.2 <- weekly.data.2$house
cue.2 <- ifelse(weekly.data.2$cue > 0, 1, 0)
inauguration <- weekly.data.2$inauguration  # 1 if an inauguration cue was provided, 0 otherwise
table(inauguration)

## estimation
# initial values
inits.inauguration <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.1)
  g <- runif(1, -0.03, 0.03)
  d <- runif(10, -0.03, 0.03)
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(1, -0.03, 0.03)
  x <- runif(1, -0.03, 0.03)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, gamma.star = g, delta.star = d, tau = t, 
       kappa = k, lambda = l, xi = x, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.inauguration <- list()
for (i in 1:3) {
  initial.values.inauguration[[i]] <- inits.inauguration(3000 + 100 * i, 4000 + 100 * i, i)
}
# MCMC
result.inauguration <- run.jags(model = "JAGS_inauguration.r", 
                                monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                            "alpha5", "alpha6", "alpha7", "alpha8", 
                                            "lambda", "xi", "gamma", "delta", "kappa", 
                                            "omega", "tau"), 
                                data = list(y = y.2, weekid = weekid.2, n.weeks = n.weeks, n.obs = n.obs.2, 
                                            survey.mode = survey.mode.2, design = design.2, house = house.2, 
                                            reshuffle = reshuffle, cue = cue.2, inauguration = inauguration, 
                                            lo = lo, up = up, start.week = start.week.2, end.week = end.week.2), 
                                inits = initial.values.inauguration, 
                                n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                                thin = 20, method = "parallel")
# save(result.inauguration, file = "result_inauguration.Rdata")
# load("result_inauguration.Rdata")

# convergence check
which(gelman.diag(result.inauguration$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.inauguration <- rbind(result.inauguration$mcmc[[1]], result.inauguration$mcmc[[2]], result.inauguration$mcmc[[3]])
colnames(sims.matrix.inauguration) <- colnames(result.inauguration$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.inauguration[, "lambda"]), 3)

## the posterior mean and the 95% HPD interval of the effect of an inauguration cue (xi)
round(posterior.summary(sims.matrix.inauguration[, "xi"]), 3)

#### Table 2 (pattern of the provision of reshuffle cues) ####
# ID of weeks in which PMs reshuffles his cabinet
reshuffle.cue.pattern <- matrix(NA, 12, 14)
for (i in 1:11) {
  # extract polls of polling firm i
  subset.weekly.data <- subset(weekly.data, house == i)
  for (j in 1:13) {
    # if polling firm i did not conduct a poll, record "-"
    if (is.na(match(reshuffle.week.id[j], subset.weekly.data$weekid) & 
              is.na(match(reshuffle.week.id[j] + 1, subset.weekly.data$weekid)))) {
      reshuffle.cue.pattern[i, j] <- "-"
    } else {
      # extract polls immediately following reshuffle j
      subset.weekly.data2 <- subset(subset.weekly.data, 
                                    weekid == reshuffle.week.id[j] | weekid == reshuffle.week.id[j] + 1)
      # if type I cue was provided
      reshuffle.cue.pattern[i, j] <- ifelse(sum(subset.weekly.data2$cue == 1) == 1, 
                                            # record "1*" if polling firm i conducted another poll, otherwise record "1"
                                            ifelse(nrow(subset.weekly.data2) == 1, "1", "1*"), 
                                            # if type II cue was provided
                                            ifelse(sum(subset.weekly.data2$cue == 2) == 1, 
                                                   # record "1_*" if polling firm i conducted another poll, otherwise record "1_"
                                                   ifelse(nrow(subset.weekly.data2) == 1, "1_", "1_*"), 
                                                   # if no cue was provided, record "0"
                                                   "0"))
    }
  }
  # total number of cues provided by polling firm i
  reshuffle.cue.pattern[i, 14] <- sum(! (reshuffle.cue.pattern[i, 1:13] == "-" | 
                                           reshuffle.cue.pattern[i, 1:13] == "0"))
}
for (j in 1:13) {
  # total number of cues provided just following reshuffle j
  reshuffle.cue.pattern[12, j] <- sum(! (reshuffle.cue.pattern[1:11, j] == "-" | 
                                           reshuffle.cue.pattern[1:11, j] == "0"))
}
# total number of cues provided by all polling firm
reshuffle.cue.pattern[12, 14] <- sum(as.numeric(reshuffle.cue.pattern[12, 1:13]))
# output Table 2
print(reshuffle.cue.pattern, quote = FALSE)

#### analysis dealing with the possibility of self-selection bias (online Appendix C) ####
## preparing the dataset, using only the data from Jiji, Kyodo, Mainichi, and NHK
weekly.data.3 <- subset(weekly.data, house == 1 | house == 2 | house == 5 | house == 8)
nrow(weekly.data.3)  # number of polls
weekid.3 <- weekly.data.3$weekid - inauguration.week[(weekly.data.3$pmid - 55)] + 1
start.week.3 <- end.week.3 <- rep(NA, 8)
counter <- 0
for (p in 1:8) {
  start.week.3[p] <- counter + 1
  end.week.3[p] <- counter <- counter + sum(weekly.data.3$pmid == (p + 55))
}
y.3 <- weekly.data.3$result / 100
n.obs.3 <- weekly.data.3$nobs
design.3 <- ifelse(weekly.data.3$mode > 2, 1, 0)
house.3 <- weekly.data.3$house
cue.3 <- ifelse(weekly.data.3$cue > 0, 1, 0)

## estimation
# initial values
inits.self.selection <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.1)
  d <- runif(7, -0.03, 0.03)
  t <- runif(8, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(1, -0.03, 0.03)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, delta.star = d, tau = t, kappa = k, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.self.selection <- list()
for (i in 1:3) {
  initial.values.self.selection[[i]] <- inits.self.selection(5000 + 100 * i, 6000 + 100 * i, i)
}
# MCMC
result.self.selection <- run.jags(model = "JAGS_self_selection.r", 
                                  monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                              "alpha5", "alpha6", "alpha7", "alpha8", 
                                              "lambda", "delta", "kappa", "omega", "tau"), 
                                  data = list(y = y.3, weekid = weekid.3, n.weeks = n.weeks, n.obs = n.obs.3, 
                                              design = design.3, house = house.3, 
                                              reshuffle = reshuffle, cue = cue.3, lo = lo, up = up, 
                                              start.week = start.week.3, end.week = end.week.3), 
                                  inits = initial.values.self.selection, 
                                  n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                                  thin = 20, method = "parallel")
# save(result.self.selection, file = "result_self_selection.Rdata")
# load("result_self_selection.Rdata")

# convergence check
which(gelman.diag(result.self.selection$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.self.selection <- rbind(result.self.selection$mcmc[[1]], result.self.selection$mcmc[[2]], result.self.selection$mcmc[[3]])
colnames(sims.matrix.self.selection) <- colnames(result.self.selection$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.self.selection[, "lambda"]), 3)

#### varying effects of reshuffle cues (online Appendix D) ####
# create indices of reshuffles
reshuffle.id <- rep(1, nrow(week.id.data))
for (i in 1:length(reshuffle.week.id)) {
  reshuffle.id[c(reshuffle.week.id[i], reshuffle.week.id[i] + 1)] <- i
}
serial.weekid <- weekly.data$weekid

# covaraites for each reshuffle
reshuffle.data <- read.csv("reshuffle_data.csv")
DPJ <- reshuffle.data$DPJ
alternation <- reshuffle.data$alternation

## estimation
# initial values
inits.varying <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.1)
  g <- runif(1, -0.03, 0.03)
  d <- runif(10, -0.03, 0.03)
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  th <- runif(4, -0.01, 0.01)
  e <- runif(1, 0.01, 0.05)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, gamma.star = g, delta.star = d, tau = t, kappa = k, 
       theta = th, eta = e, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.varying <- list()
for (i in 1:3) {
  initial.values.varying[[i]] <- inits.varying(7000 + 100 * i, 8000 + 100 * i, i)
}
# MCMC
result.varying <- run.jags(model = "JAGS_varying_effects.r", 
                           monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                       "alpha5", "alpha6", "alpha7", "alpha8", 
                                       "lambda", "theta", "gamma", "delta", "kappa", 
                                       "omega", "tau", "eta", "approval"), 
                           data = list(y = y, weekid = weekid, serial.weekid = serial.weekid, 
                                       n.weeks = n.weeks, n.obs = n.obs, 
                                       survey.mode = survey.mode, design = design, house = house, 
                                       reshuffle = reshuffle, cue = cue, lo = lo, up = up, 
                                       start.week = start.week, end.week = end.week, 
                                       reshuffle.id = reshuffle.id, DPJ = DPJ, alternation = alternation), 
                           inits = initial.values.varying, 
                           n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                           thin = 20, method = "parallel")
# save(result.varying, file = "result_varying.Rdata")
# load("result_varying.Rdata")

# convergence check
which(gelman.diag(result.varying$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.varying <- rbind(result.varying$mcmc[[1]], result.varying$mcmc[[2]], result.varying$mcmc[[3]])
colnames(sims.matrix.varying) <- colnames(result.varying$mcmc[[1]])

## Table 3 (the posterior mean and the 95% HPD interval of coefficients)
varying.theta.result <- matrix(NA, 4, 3)
for (i in 1:4) {
  varying.theta.result[i, ] <- posterior.summary(sims.matrix.varying[, 775 + i])
}
round(varying.theta.result, 3)

# 80% CI of the coefficient of the number of altered ministers does not include zero
0 < HPDinterval(as.mcmc(sims.matrix.varying[, "theta[3]"]), prob = 0.8)[1] | 
  HPDinterval(as.mcmc(sims.matrix.varying[, "theta[3]"]), prob = 0.8)[2] < 0
round(HPDinterval(as.mcmc(sims.matrix.varying[, "theta[3]"]), prob = 0.8), 3)

## Table S2
round(cbind(colMeans(sims.matrix.varying[, 814:826]), alternation, DPJ), 3)

## Figure S1 (predicted values of the size of the effect of a reshuffle cue for each reshuffle)
varying.lambda.result <- matrix(NA, 13, 3)
for (i in 1:13) {
  varying.lambda.result[i, ] <- posterior.summary(sims.matrix.varying[, 762 + i])
}
round(varying.lambda.result, 3)

png("Figure_S1.png", width = 120, height = 100, units = "mm", pointsize = 8, res = 800)
par(mar = c(5, 2, 2, 2))
dotchart(rev(varying.lambda.result[, 1]), 
         labels = c("Abe (2), Oct 2015", "Abe (2), Sep 2014", "Noda, Oct 2012", "Noda, Jun 2012", 
                    "Noda, Jan 2012", "Kan, Jan 2011", "Kan, Sep 2010", "Fukuda, Aug 2008", 
                    "Abe (1), Aug 2007", "Koizumi, Oct 2005", "Koizumi, Sep 2004", 
                    "Koizumi, Sep 2003", "Koizumi, Sep 2002"), 
         pch = 19, xlim = c(-0.1, 0.1), xlab = "Effects (Point)", ylab = "")
points(varying.lambda.result[, 1], c(13:1), pch = 19, col = color[12])
segments(varying.lambda.result[, 2], c(13:1), varying.lambda.result[, 3], c(13:1), col = color[12])
abline(v = 0, lty = 3)
dev.off()

#### analysis of the two types of reshuffle cues (online Appendix E) ####
## variables for two types of cue
cue.1 <- ifelse(weekly.data$cue == 1, 1, 0)  # type I
cue.2 <- ifelse(weekly.data$cue == 2, 1, 0)  # type II

## estimation
# initial values
inits.two.types <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.1)
  g <- runif(1, -0.03, 0.03)
  d <- runif(10, -0.03, 0.03)
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(2, -0.03, 0.03)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, gamma.star = g, delta.star = d, tau = t, kappa = k, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.two.types <- list()
for (i in 1:3) {
  initial.values.two.types[[i]] <- inits.two.types(9000 + 100 * i, 10000 + 100 * i, i)
}
# MCMC
result.two.types <- run.jags(model = "JAGS_two_types.r", 
                             monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                         "alpha5", "alpha6", "alpha7", "alpha8", 
                                         "lambda", "gamma", "delta", "kappa", 
                                         "omega", "tau"), 
                             data = list(y = y, weekid = weekid, n.weeks = n.weeks, n.obs = n.obs, 
                                         survey.mode = survey.mode, design = design, house = house, 
                                         reshuffle = reshuffle, cue.1 = cue.1, cue.2 = cue.2, lo = lo, up = up, 
                                         start.week = start.week, end.week = end.week), 
                             inits = initial.values.two.types, 
                             n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                             thin = 20, method = "parallel")
# save(result.two.types, file = "result_two_types.Rdata")
# load("result_two_types.Rdata")

# convergence check
which(gelman.diag(result.two.types$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.two.types <- rbind(result.two.types$mcmc[[1]], result.two.types$mcmc[[2]], result.two.types$mcmc[[3]])
colnames(sims.matrix.two.types) <- colnames(result.two.types$mcmc[[1]])

## Figure S2 (estimated effect of a reshuffle cue according to type of cues)
# extracting parameter samples
lambda1.samples <- sims.matrix.two.types[, "lambda[1]"]
lambda2.samples <- sims.matrix.two.types[, "lambda[2]"]
# the posterior means and the 95% HPD intervals
post.lambda1.mean <- mean(lambda1.samples)
post.lambda1.HPD <- HPDinterval(as.mcmc(lambda1.samples))
post.lambda2.mean <- mean(lambda2.samples)
post.lambda2.HPD <- HPDinterval(as.mcmc(lambda2.samples))
# difference in two types
lambda.dif.samples <- sims.matrix.two.types[, "lambda[2]"] - sims.matrix.two.types[, "lambda[1]"]
round(posterior.summary(lambda.dif.samples), 3)

# drawing Figure S2
png("Figure_S2.png", width = 120, height = 50, units = "mm", pointsize = 8, res = 800)
par(mar = c(5, 2, 2, 2))
dotchart(c(post.lambda2.mean, post.lambda1.mean, post.lambda.mean), 
         labels = c("Type II", "Type I", "Combined"), 
         pch = 19, xlim = c(0, 0.06), xlab = "Effects (Point)", ylab = "")
points(post.lambda.mean, 3, pch = 19, col = color[12])
points(post.lambda1.mean, 2, pch = 19, col = color[2])
points(post.lambda2.mean, 1, pch = 19, col = color[2])
segments(post.lambda.HPD[1], 3, post.lambda.HPD[2], 3, col = color[12])
segments(post.lambda1.HPD[1], 2, post.lambda1.HPD[2], 2, col = color[2])
segments(post.lambda2.HPD[1], 1, post.lambda2.HPD[2], 1, col = color[2])
abline(v = 0, lty = 3)
dev.off()

#### robustness check with a simpler model (online Appendix F) ####
## preparing data of polls conducted just after reshuffles
reshuffle.week.data <- subset(weekly.data, reshuffle == 1)
reshuffle.week.data$result <- reshuffle.week.data$result / 100
reshuffle.episode.label <- c("01.Koizumi_1", "02.Koizumi_2", "03.Koizumi_3", "04.Koizumi_4", 
                             "05.Abe_1", "06.Fukuda", "07.Kan_1", "08.Kan_2", "09.Noda_1", "10.Noda_2", 
                             "11.Noda_3", "12.Abe_2", "13.Abe_3")
reshuffle.episode <- rep(NA, nrow(reshuffle.week.data))  # factors of each reshuffle
for (i in 1:nrow(reshuffle.week.data)) {
  for (j in 1:length(reshuffle.episode.label)) {
    if (reshuffle.week.data$weekid[i] == reshuffle.week.id[j] | reshuffle.week.data$weekid[i] == reshuffle.week.id[j] + 1) {
      reshuffle.episode[i] <- reshuffle.episode.label[j]
      break
    }
  }
}
reshuffle.week.data <- data.frame(reshuffle.week.data, reshuffle.episode)

## estimating by OLS without any corrections for an estimated dependent variable
lm.result <- lm(result ~ reshuffle.episode + factor(house) + I(cue > 0) - 1, 
                data = reshuffle.week.data, x = TRUE)

## OLS with Efron�fs SE estimator (first column of Table S3)
round(coeftest(lm.result, vcov. = vcovHC, type = "HC3"), 3)

## FGLS assuming that the sampling error is known (second column of Table S3)
sq.residuals <- lm.result$residuals ^ 2
sq.omega <- reshuffle.week.data$result * (1 - reshuffle.week.data$result) / reshuffle.week.data$nobs
G <- diag(reshuffle.week.data$nobs)
X <- lm.result$x
sq.sigma.hat.1 <- (sum(sq.residuals) - sum(sq.omega) + 
                     sum(diag(solve(crossprod(X) %*% t(X) %*% G %*% X)))) / lm.result$df.residual
fglm.weights.1 <- 1 / sqrt(sq.omega + sq.sigma.hat.1)
fglm.result.1 <- lm(result ~ reshuffle.episode + factor(house) + I(cue > 0) - 1, 
                    data = reshuffle.week.data, weights = fglm.weights.1)
round(coeftest(fglm.result.1), 3)

## FGLS assuming that the sampling error is proportional (third column of Table S3)
lm.kappa <- lm(sq.residuals ~ sq.omega)
sq.sigma.hat.2 <- lm.kappa$coefficients[1]
kappa.hat <- lm.kappa$coefficients[2]
fglm.weights.2 <- 1 / sqrt(sq.sigma.hat.2 + kappa.hat * sq.omega)
fglm.result.2 <- lm(result ~ reshuffle.episode + factor(house) + I(cue > 0) - 1, 
                    data = reshuffle.week.data, weights = fglm.weights.2)
round(coeftest(fglm.result.2), 3)

#### AR(1) model ####
## estimation
# initial values
inits.AR1 <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.1)
  g <- runif(1, -0.03, 0.03)
  d <- runif(10, -0.03, 0.03)
  ph <- runif(8, 0.9, 1)
  ps <- 0.4 * (1 - ph)
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(1, -0.03, 0.03)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, gamma.star = g, delta.star = d, psi = ps, phi = ph, tau = t, kappa = k, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.AR1 <- list()
for (i in 1:3) {
  initial.values.AR1[[i]] <- inits.AR1(11000 + 100 * i, 12000 + 100 * i, i)
}
# MCMC
result.AR1 <- run.jags(model = "JAGS_AR1.r", 
                       monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                   "alpha5", "alpha6", "alpha7", "alpha8", 
                                   "lambda", "gamma", "delta", "kappa", 
                                   "psi", "phi", "omega", "tau"), 
                       data = list(y = y, weekid = weekid, n.weeks = n.weeks, n.obs = n.obs, 
                                   survey.mode = survey.mode, design = design, house = house, 
                                   reshuffle = reshuffle, cue = cue, lo = lo, up = up, 
                                   start.week = start.week, end.week = end.week), 
                       inits = initial.values.AR1, 
                       n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                       thin = 20, method = "parallel")
# save(result.AR1, file = "result_AR1.Rdata")
# load("result_AR1.Rdata")

# convergence check
which(gelman.diag(result.AR1$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.AR1 <- rbind(result.AR1$mcmc[[1]], result.AR1$mcmc[[2]], result.AR1$mcmc[[3]])
colnames(sims.matrix.AR1) <- colnames(result.AR1$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.AR1[, "lambda"]), 3)

#### AR(2) Model ####
## estimation
# initial values
inits.AR2 <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.1)
  g <- runif(1, -0.03, 0.03)
  d <- runif(10, -0.03, 0.03)
  ph <- runif(16, 0.9, 1)
  ps <- 0.4 * (1 - ph[1:8] - ph[9:16])
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(1, -0.03, 0.03)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, gamma.star = g, delta.star = d, psi = ps, phi = matrix(ph, 8, 2), tau = t, kappa = k, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.AR2 <- list()
for (i in 1:3) {
  initial.values.AR2[[i]] <- inits.AR2(13000 + 100 * i, 14000 + 100 * i, i)
}
# MCMC
result.AR2 <- run.jags(model = "JAGS_AR2.r", 
                       monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                   "alpha5", "alpha6", "alpha7", "alpha8", 
                                   "lambda", "gamma", "delta", "kappa", 
                                   "psi", "phi", "omega", "tau"), 
                       data = list(y = y, weekid = weekid, n.weeks = n.weeks, n.obs = n.obs, 
                                   survey.mode = survey.mode, design = design, house = house, 
                                   reshuffle = reshuffle, cue = cue, lo = lo, up = up, 
                                   start.week = start.week, end.week = end.week), 
                       inits = initial.values.AR2, 
                       n.chains = 3, burnin = 5000, sample = 2000, adapt = 5000, summarise = FALSE, 
                       thin = 150, method = "parallel")
# save(result.AR2, file = "result_AR2.Rdata")
# load("result_AR2.Rdata")

# convergence check
which(gelman.diag(result.AR2$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.AR2 <- rbind(result.AR2$mcmc[[1]], result.AR2$mcmc[[2]], result.AR2$mcmc[[3]])
colnames(sims.matrix.AR2) <- colnames(result.AR2$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.AR2[, "lambda"]), 3)

#### analysis without Kyodo News data ####
## preparing dataset
weekly.data.4 <- subset(weekly.data, house != 2)  # excluding Kyodo News data
nrow(weekly.data.4)  # number of polls
weekid.4 <- weekly.data.4$weekid - inauguration.week[(weekly.data.4$pmid - 55)] + 1
start.week.4 <- end.week.4 <- rep(NA, 8)
counter <- 0
for (p in 1:8) {
  start.week.4[p] <- counter + 1
  end.week.4[p] <- counter <- counter + sum(weekly.data.4$pmid == (p + 55))
}
y.4 <- weekly.data.4$result / 100
n.obs.4 <- weekly.data.4$nobs
survey.mode.4 <- ifelse(weekly.data.4$mode == 1, 1, 2)
design.4 <- ifelse(weekly.data.4$mode > 2, 1, 0)
house.4 <- weekly.data.4$house
cue.4 <- ifelse(weekly.data.4$cue > 0, 1, 0)

## estimation
# preparing three sets of initial values for three MCMC chains
initial.values.wo.Kyodo <- list()
for (i in 1:3) {
  initial.values.wo.Kyodo[[i]] <- inits.main(15000 + 100 * i, 16000 + 100 * i, i)
}
# MCMC
result.wo.Kyodo <- run.jags(model = "JAGS_wo_Kyodo.r", 
                            monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                        "alpha5", "alpha6", "alpha7", "alpha8", 
                                        "lambda", "gamma", "delta", "kappa", 
                                        "psi", "phi", "omega", "tau"), 
                            data = list(y = y.4, weekid = weekid.4, n.weeks = n.weeks, n.obs = n.obs.4, 
                                        survey.mode = survey.mode.4, design = design.4, house = house.4, 
                                        reshuffle = reshuffle, cue = cue.4, lo = lo, up = up, 
                                        start.week = start.week.4, end.week = end.week.4), 
                            inits = initial.values.wo.Kyodo, 
                            n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                            thin = 20, method = "parallel")
# save(result.wo.Kyodo, file = "result_wo_Kyodo.Rdata")
# load("result_wo_Kyodo.Rdata")

# convergence check
which(gelman.diag(result.wo.Kyodo$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.wo.Kyodo <- rbind(result.wo.Kyodo$mcmc[[1]], result.wo.Kyodo$mcmc[[2]], result.wo.Kyodo$mcmc[[3]])
colnames(sims.matrix.wo.Kyodo) <- colnames(result.wo.Kyodo$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.wo.Kyodo[, "lambda"]), 3)

#### analysis of whether house effects vary by government partisanship ####
## preparing the variable for government partisanship
# party = 1 if the government is the LDP's coalition
# party = 2 if the government is the DPJ's coalition
party <- ifelse(weekly.data$pmid == 60 | weekly.data$pmid == 61 | weekly.data$pmid == 62, 2, 1)

## estimation
# initial values
inits.party <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.1)
  g <- runif(1, -0.03, 0.03)
  d <- runif(20, -0.03, 0.03)
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(1, -0.03, 0.03)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, gamma.star = g, delta.star = matrix(d, 10, 2), 
       tau = t, kappa = k, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.party <- list()
for (i in 1:3) {
  initial.values.party[[i]] <- inits.party(17000 + 100 * i, 18000 + 100 * i, i)
}
# MCMC
result.party <- run.jags(model = "JAGS_party.r", 
                         monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                     "alpha5", "alpha6", "alpha7", "alpha8", 
                                     "lambda", "gamma", "delta", "kappa", 
                                     "omega", "tau"), 
                         data = list(y = y, weekid = weekid, n.weeks = n.weeks, n.obs = n.obs, 
                                     survey.mode = survey.mode, design = design, house = house, 
                                     reshuffle = reshuffle, cue = cue, party = party, lo = lo, up = up, 
                                     start.week = start.week, end.week = end.week), 
                         inits = initial.values.party, 
                         n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                         thin = 20, method = "parallel")
# save(result.party, file = "result_party.Rdata")
# load("result_party.Rdata")

# convergence check
which(gelman.diag(result.party$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.party <- rbind(result.party$mcmc[[1]], result.party$mcmc[[2]], result.party$mcmc[[3]])
colnames(sims.matrix.party) <- colnames(result.party$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.party[, "lambda"]), 3)

## Figure S3 (estimated effect of house effects depending on government partisanship)
# extracting parameters
delta.samples.LDP <- sims.matrix.party[, 766:776]
delta.samples.DPJ <- sims.matrix.party[, 777:787]

# computing the posterior means and 95% HPD intervals
post.delta.LDP.mean <- colMeans(delta.samples.LDP)
post.delta.LDP.HPD <- HPDinterval(as.mcmc(delta.samples.LDP))

post.delta.DPJ.mean <- colMeans(delta.samples.DPJ)
post.delta.DPJ.HPD <- HPDinterval(as.mcmc(delta.samples.DPJ))

# drawing Figure S3
png("Figure_S3.png", width = 120, height = 100, units = "mm", pointsize = 8, res = 800)
par(mar = c(5, 2, 2, 2))
dotchart(rep(1, 11), 
         labels = rev(c("Jiji", "Kyodo", "Yomiuri", "Asahi", "Mainichi", 
                        "Sankei-FNN", "Nikkei", "NHK", "JNN", "ANN", "NNN")), 
         pch = 19, xlim = c(-0.07, 0.07), xlab = "Effects (Point)", ylab = "")
points(post.delta.LDP.mean, c(11:1) + 0.05, pch = 19, col = color[12])
points(post.delta.DPJ.mean, c(11:1) - 0.05, pch = 15, col = color[10])
segments(post.delta.LDP.HPD[, 1], c(11:1) + 0.05, post.delta.LDP.HPD[, 2], c(11:1) + 0.05, col = color[12])
segments(post.delta.DPJ.HPD[, 1], c(11:1) - 0.05, post.delta.DPJ.HPD[, 2], c(11:1) - 0.05, col = color[10])
abline(v = 0, lty = 3)
dev.off()

# polling firms whose house effect vary considerably depending on government partisanship
round(mean(delta.samples.DPJ[, 5] - delta.samples.LDP[, 5]), 3)  # Mainichi
round(mean(delta.samples.DPJ[, 7] - delta.samples.LDP[, 7]), 3)  # Nikkei
round(mean(delta.samples.DPJ[, 9] - delta.samples.LDP[, 9]), 3)  # JNN

#### analysis altering the coding of Nikkei Research�fs poll conducted in September 2010 ####
## preparing variables
cue.Nikkei <- cue
# recoding from a Type I cue (cue = 1) to no cues (cue = 0)
cue.Nikkei[house == 7 & weekid == 491] <- 0

## estimation
# preparing three sets of initial values for three MCMC chains
initial.values.Nikkei1 <- list()
for (i in 1:3) {
  initial.values.Nikkei1[[i]] <- inits.main(19000 + 100 * i, 20000 + 100 * i, i)
}
# MCMC
result.Nikkei1 <- run.jags(model = "JAGS_main.r", 
                       monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                   "alpha5", "alpha6", "alpha7", "alpha8", 
                                   "lambda", "gamma", "delta", "kappa", 
                                   "omega", "tau"), 
                       data = list(y = y, weekid = weekid, n.weeks = n.weeks, n.obs = n.obs, 
                                   survey.mode = survey.mode, design = design, house = house, 
                                   reshuffle = reshuffle, cue = cue.Nikkei, lo = lo, up = up, 
                                   start.week = start.week, end.week = end.week), 
                       inits = initial.values.Nikkei1, 
                       n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                       thin = 20, method = "parallel")
# save(result.Nikkei1, file = "result_Nikkei1.Rdata")
# load("result_Nikkei1.Rdata")

# convergence check
which(gelman.diag(result.Nikkei1$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.Nikkei1 <- rbind(result.Nikkei1$mcmc[[1]], result.Nikkei1$mcmc[[2]], result.Nikkei1$mcmc[[3]])
colnames(sims.matrix.Nikkei1) <- colnames(result.Nikkei1$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.Nikkei1[, "lambda"]), 3)

#### analysis omitting Nikkei Research�fs poll conducted in September 2010 ####
## preparing dataset
weekly.data.5 <- subset(weekly.data, !(house == 7 & weekid == 491))  # excludin the concerned poll
nrow(weekly.data.5)  # number of polls
weekid.5 <- weekly.data.5$weekid - inauguration.week[(weekly.data.5$pmid - 55)] + 1
start.week.5 <- end.week.5 <- rep(NA, 8)
counter <- 0
for (p in 1:8) {
  start.week.5[p] <- counter + 1
  end.week.5[p] <- counter <- counter + sum(weekly.data.5$pmid == (p + 55))
}
y.5 <- weekly.data.5$result / 100
n.obs.5 <- weekly.data.5$nobs
survey.mode.5 <- ifelse(weekly.data.5$mode == 1, 1, 2)
design.5 <- ifelse(weekly.data.5$mode > 2, 1, 0)
house.5 <- weekly.data.5$house
cue.5 <- ifelse(weekly.data.5$cue > 0, 1, 0)

## estimation
# preparing three sets of initial values for three MCMC chains
initial.values.Nikkei2 <- list()
for (i in 1:3) {
  initial.values.Nikkei2[[i]] <- inits.main(21000 + 100 * i, 22000 + 100 * i, i)
}
# MCMC
result.Nikkei2 <- run.jags(model = "JAGS_main.r", 
                           monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                       "alpha5", "alpha6", "alpha7", "alpha8", 
                                       "lambda", "gamma", "delta", "kappa", 
                                       "omega", "tau"), 
                           data = list(y = y.5, weekid = weekid.5, n.weeks = n.weeks, n.obs = n.obs.5, 
                                       survey.mode = survey.mode.5, design = design.5, house = house.5, 
                                       reshuffle = reshuffle, cue = cue.5, lo = lo, up = up, 
                                       start.week = start.week.5, end.week = end.week.5), 
                           inits = initial.values.Nikkei2, 
                           n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                           thin = 20, method = "parallel")
# save(result.Nikkei2, file = "result_Nikkei2.Rdata")
# load("result_Nikkei2.Rdata")

# convergence check
which(gelman.diag(result.Nikkei2$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.Nikkei2 <- rbind(result.Nikkei2$mcmc[[1]], result.Nikkei2$mcmc[[2]], result.Nikkei2$mcmc[[3]])
colnames(sims.matrix.Nikkei2) <- colnames(result.Nikkei2$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.Nikkei2[, "lambda"]), 3)

#### analysis altering the mode of NNN�fs opinion polls before December 2006 ####
## preparing variables
survey.mode.NNN <- survey.mode
# recoding from telephone interviews (survey.mode = 2) to face-to-face interview (survey.mode = 1)
survey.mode.NNN[house == 11 & survey.mode == 2 & design == 0] <- 1

## estimation
# preparing three sets of initial values for three MCMC chains
initial.values.NNN <- list()
for (i in 1:3) {
  initial.values.NNN[[i]] <- inits.main(23000 + 100 * i, 24000 + 100 * i, i)
}
# MCMC estimation
result.NNN <- run.jags(model = "JAGS_main.r", 
                       monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                   "alpha5", "alpha6", "alpha7", "alpha8", 
                                   "lambda", "gamma", "delta", "kappa", 
                                   "omega", "tau"), 
                       data = list(y = y, weekid = weekid, n.weeks = n.weeks, n.obs = n.obs, 
                                   survey.mode = survey.mode.NNN, design = design, house = house, 
                                   reshuffle = reshuffle, cue = cue, lo = lo, up = up, 
                                   start.week = start.week, end.week = end.week), 
                       inits = initial.values.NNN, 
                       n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                       thin = 20, method = "parallel")
# save(result.NNN, file = "result_NNN.Rdata")
# load("result_NNN.Rdata")

# convergence check
which(gelman.diag(result.NNN$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.NNN <- rbind(result.NNN$mcmc[[1]], result.NNN$mcmc[[2]], result.NNN$mcmc[[3]])
colnames(sims.matrix.NNN) <- colnames(result.NNN$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.NNN[, "lambda"]), 3)

#### analysis without the reshuffle dummy ####
## estimation
# initial values
inits.wo.kappa <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- rnorm(n.weeks[1], lo[1] + 0.1, 0.01)
  a2 <- rnorm(n.weeks[2], lo[2] + 0.1, 0.01)
  a3 <- rnorm(n.weeks[3], lo[3] + 0.1, 0.01)
  a4 <- rnorm(n.weeks[4], lo[4] + 0.1, 0.01)
  a5 <- rnorm(n.weeks[5], lo[5] + 0.1, 0.01)
  a6 <- rnorm(n.weeks[6], lo[6] + 0.1, 0.01)
  a7 <- rnorm(n.weeks[7], lo[7] + 0.1, 0.01)
  a8 <- rnorm(n.weeks[8], lo[8] + 0.1, 0.01)
  o <- runif(8, 0.01, 0.5)
  g <- runif(1, -0.03, 0.03)
  d <- runif(10, -0.03, 0.03)
  t <- runif(11, 1, 2)
  l <- runif(1, -0.03, 0.03)
  list(alpha1 = a1, alpha2 = a2, alpha3 = a3, alpha4 = a4, alpha5 = a5, 
       alpha6 = a6, alpha7 = a7, alpha8 = a8, 
       omega = o, gamma.star = g, delta.star = d, tau = t, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.wo.kappa <- list()
for (i in 1:3) {
  initial.values.wo.kappa[[i]] <- inits.wo.kappa(25000 + 100 * i, 26000 + 100 * i, i)
}
# MCMC
result.wo.kappa <- run.jags(model = "JAGS_wo_kappa.r", 
                            monitor = c("alpha1", "alpha2", "alpha3", "alpha4", 
                                        "alpha5", "alpha6", "alpha7", "alpha8", 
                                        "lambda", "gamma", "delta", "omega", "tau"), 
                            data = list(y = y, weekid = weekid, n.weeks = n.weeks, n.obs = n.obs, 
                                        survey.mode = survey.mode, design = design, house = house, 
                                        cue = cue, lo = lo, up = up, 
                                        start.week = start.week, end.week = end.week), 
                            inits = initial.values.wo.kappa, 
                            n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                            thin = 20, method = "parallel")
# save(result.wo.kappa, file = "result_wo_kappa.Rdata")
# load("result_wo_kappa.Rdata")

# convergence check
which(gelman.diag(result.wo.kappa$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.wo.kappa <- rbind(result.wo.kappa$mcmc[[1]], result.wo.kappa$mcmc[[2]], result.wo.kappa$mcmc[[3]])
colnames(sims.matrix.wo.kappa) <- colnames(result.wo.kappa$mcmc[[1]])

## the posterior mean and the 95% HPD interval of the effect of a reshuffle cue (lambda)
round(posterior.summary(sims.matrix.wo.kappa[, "lambda"]), 3)

#### Analysis for each PM (online Appendix G) ####
## Koizumi
# initial values
inits.single.PM <- function (n.weeks, lo, seed1, seed2, i) {
  set.seed(seed1)
  a <- rnorm(n.weeks, lo + 0.1, 0.01)
  o <- runif(1, 0.01, 0.1)
  d <- runif(10, -0.03, 0.03)
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(1, -0.03, 0.03)
  list(alpha = a, omega = o, delta.star = d, tau = t, kappa = k, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.Koizumi <- list()
for (i in 1:3) {
  initial.values.Koizumi[[i]] <- inits.single.PM(n.weeks = n.weeks[1], lo = lo[1], 
                                                 27000 + 100 * i, 28000 + 100 * i, i)
}
# MCMC
result.Koizumi <- run.jags(model = "JAGS_Koizumi.r", 
                           monitor = c("alpha", "lambda", "delta", "kappa", "omega", "tau"), 
                           data = list(y = y[weekly.data$pmid == 56], weekid = weekid[weekly.data$pmid == 56], 
                                       n.weeks = n.weeks[1], n.obs = n.obs[weekly.data$pmid == 56], 
                                       design = design[weekly.data$pmid == 56], 
                                       house = house[weekly.data$pmid == 56], 
                                       reshuffle = reshuffle[1, ], 
                                       cue = cue[weekly.data$pmid == 56], 
                                       lo = lo[1], up = up[1], n.polls = length(y[weekly.data$pmid == 56])), 
                           inits = initial.values.Koizumi, 
                           n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                           thin = 20, method = "parallel")
# save(result.Koizumi, file = "result_Koizumi.Rdata")
# load("result_Koizumi.Rdata")

# convergence check
which(gelman.diag(result.Koizumi$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.Koizumi <- rbind(result.Koizumi$mcmc[[1]], result.Koizumi$mcmc[[2]], result.Koizumi$mcmc[[3]])
colnames(sims.matrix.Koizumi) <- colnames(result.Koizumi$mcmc[[1]])

## Abe (1)
# preparing three sets of initial values for three MCMC chains
initial.values.Abe1 <- list()
for (i in 1:3) {
  initial.values.Abe1[[i]] <- inits.single.PM(n.weeks = n.weeks[2], lo = lo[2], 
                                              29000 + 100 * i, 30000 + 100 * i, i)
}
# MCMC
result.Abe1 <- run.jags(model = "JAGS_each_PM.r", 
                        monitor = c("alpha", "lambda", "delta", "kappa", "omega", "tau"), 
                        data = list(y = y[weekly.data$pmid == 57], weekid = weekid[weekly.data$pmid == 57], 
                                    n.weeks = n.weeks[2], n.obs = n.obs[weekly.data$pmid == 57], 
                                    design = design[weekly.data$pmid == 57], 
                                    house = house[weekly.data$pmid == 57], 
                                    reshuffle = reshuffle[2, ], 
                                    cue = cue[weekly.data$pmid == 57], 
                                    lo = lo[2], up = up[2], n.polls = length(y[weekly.data$pmid == 57])), 
                        inits = initial.values.Abe1, 
                        n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                        thin = 20, method = "parallel")
# save(result.Abe1, file = "result_Abe1.Rdata")
# load("result_Abe1.Rdata")

# convergence check
which(gelman.diag(result.Abe1$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.Abe1 <- rbind(result.Abe1$mcmc[[1]], result.Abe1$mcmc[[2]], result.Abe1$mcmc[[3]])
colnames(sims.matrix.Abe1) <- colnames(result.Abe1$mcmc[[1]])

## Fukuda
# preparing three sets of initial values for three MCMC chains
initial.values.Fukuda <- list()
for (i in 1:3) {
  initial.values.Fukuda[[i]] <- inits.single.PM(n.weeks = n.weeks[3], lo = lo[3], 
                                                31000 + 100 * i, 32000 + 100 * i, i)
}
# MCMC
result.Fukuda <- run.jags(model = "JAGS_each_PM.r", 
                          monitor = c("alpha", "lambda", "delta", "kappa", "omega", "tau"), 
                          data = list(y = y[weekly.data$pmid == 58], weekid = weekid[weekly.data$pmid == 58], 
                                      n.weeks = n.weeks[3], n.obs = n.obs[weekly.data$pmid == 58],
                                      design = design[weekly.data$pmid == 58], 
                                      house = house[weekly.data$pmid == 58], 
                                      reshuffle = reshuffle[3, ], 
                                      cue = cue[weekly.data$pmid == 58], 
                                      lo = lo[3], up = up[3], n.polls = length(y[weekly.data$pmid == 58])), 
                          inits = initial.values.Fukuda, 
                          n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                          thin = 20, method = "parallel")
# save(result.Fukuda, file = "result_Fukuda.Rdata")
# load("result_Fukuda.Rdata")

# convergence check
which(gelman.diag(result.Fukuda$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.Fukuda <- rbind(result.Fukuda$mcmc[[1]], result.Fukuda$mcmc[[2]], result.Fukuda$mcmc[[3]])
colnames(sims.matrix.Fukuda) <- colnames(result.Fukuda$mcmc[[1]])

## Kan
# initial values
inits.single.PM2 <- function (n.weeks, lo, seed1, seed2, i) {
  set.seed(seed1)
  a <- rnorm(n.weeks, lo + 0.1, 0.01)
  o <- runif(1, 0.01, 0.1)
  d <- runif(10, -0.03, 0.03)
  t <- runif(11, 1, 2)
  k <- runif(1, -0.03, 0.03)
  l <- runif(1, -0.03, 0.03)
  list(alpha = a, omega = o, delta.star = d, tau = t, kappa = k, lambda = l, 
       .RNG.name = ifelse(i == 1, "base::Wichmann-Hill", 
                          ifelse(i == 2, "base::Marsaglia-Multicarry", "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values.Kan <- list()
for (i in 1:3) {
  initial.values.Kan[[i]] <- inits.single.PM2(n.weeks = n.weeks[6], lo = lo[6], 
                                              33000 + 100 * i, 34000 + 100 * i, i)
}
# MCMC
result.Kan <- run.jags(model = "JAGS_each_PM.r", 
                       monitor = c("alpha", "lambda", "delta", "kappa", "omega", "tau"), 
                       data = list(y = y[weekly.data$pmid == 61], weekid = weekid[weekly.data$pmid == 61], 
                                   n.weeks = n.weeks[6], n.obs = n.obs[weekly.data$pmid == 61], 
                                   design = design[weekly.data$pmid == 61], 
                                   house = house[weekly.data$pmid == 61], 
                                   reshuffle = reshuffle[6, ], 
                                   cue = cue[weekly.data$pmid == 61], 
                                   lo = lo[6], up = up[6], n.polls = length(y[weekly.data$pmid == 61])), 
                       inits = initial.values.Kan, 
                       n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                       thin = 20, method = "parallel")
# save(result.Kan, file = "result_Kan.Rdata")
# load("result_Kan.Rdata")

# convergence check
which(gelman.diag(result.Kan$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.Kan <- rbind(result.Kan$mcmc[[1]], result.Kan$mcmc[[2]], result.Kan$mcmc[[3]])
colnames(sims.matrix.Kan) <- colnames(result.Kan$mcmc[[1]])

## Noda
# preparing three sets of initial values for three MCMC chains
initial.values.Noda <- list()
for (i in 1:3) {
  initial.values.Noda[[i]] <- inits.single.PM2(n.weeks = n.weeks[7], lo = lo[7], 
                                               35000 + 100 * i, 36000 + 100 * i, i)
}
# MCMC
result.Noda <- run.jags(model = "JAGS_each_PM.r", 
                        monitor = c("alpha", "lambda", "delta", "kappa", "omega", "tau"), 
                        data = list(y = y[weekly.data$pmid == 62], weekid = weekid[weekly.data$pmid == 62], 
                                    n.weeks = n.weeks[7], n.obs = n.obs[weekly.data$pmid == 62], 
                                    design = design[weekly.data$pmid == 62], 
                                    house = house[weekly.data$pmid == 62], 
                                    reshuffle = reshuffle[7, ], 
                                    cue = cue[weekly.data$pmid == 62], 
                                    lo = lo[7], up = up[7], n.polls = length(y[weekly.data$pmid == 62])), 
                        inits = initial.values.Noda, 
                        n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                        thin = 20, method = "parallel")
# save(result.Noda, file = "result_Noda.Rdata")
# load("result_Noda.Rdata")

# convergence check
which(gelman.diag(result.Noda$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.Noda <- rbind(result.Noda$mcmc[[1]], result.Noda$mcmc[[2]], result.Noda$mcmc[[3]])
colnames(sims.matrix.Noda) <- colnames(result.Noda$mcmc[[1]])

## Abe2
# preparing three sets of initial values for three MCMC chains
initial.values.Abe2 <- list()
for (i in 1:3) {
  initial.values.Abe2[[i]] <- inits.single.PM2(n.weeks = n.weeks[8], lo = lo[8], 
                                               37000 + 100 * i, 38000 + 100 * i, i)
}
# MCMC
result.Abe2 <- run.jags(model = "JAGS_each_PM.r", 
                        monitor = c("alpha", "lambda", "delta", "kappa", "omega", "tau"), 
                        data = list(y = y[weekly.data$pmid == 63], weekid = weekid[weekly.data$pmid == 63], 
                                    n.weeks = n.weeks[8], n.obs = n.obs[weekly.data$pmid == 63], 
                                    design = design[weekly.data$pmid == 63], 
                                    house = house[weekly.data$pmid == 63], 
                                    reshuffle = reshuffle[8, ], 
                                    cue = cue[weekly.data$pmid == 63], 
                                    lo = lo[8], up = up[8], n.polls = length(y[weekly.data$pmid == 63])), 
                        inits = initial.values.Abe2, 
                        n.chains = 3, burnin = 500, sample = 2000, adapt = 500, summarise = FALSE, 
                        thin = 20, method = "parallel")
# save(result.Abe2, file = "result_Abe2.Rdata")
# load("result_Abe2.Rdata")

# convergence check
which(gelman.diag(result.Abe2$mcmc, multivariate = FALSE)[[1]][, 1] > 1.1)
# conbining three chains into a single matrix
sims.matrix.Abe2 <- rbind(result.Abe2$mcmc[[1]], result.Abe2$mcmc[[2]], result.Abe2$mcmc[[3]])
colnames(sims.matrix.Abe2) <- colnames(result.Abe2$mcmc[[1]])

## Figure S4 (effect of a reshuffle cue for each PM)
each.PM.lambda.result <- matrix(NA, 6, 3)
each.PM.lambda.result[1, ] <- posterior.summary(sims.matrix.Koizumi[, "lambda"])
each.PM.lambda.result[2, ] <- posterior.summary(sims.matrix.Abe1[, "lambda"])
each.PM.lambda.result[3, ] <- posterior.summary(sims.matrix.Fukuda[, "lambda"])
each.PM.lambda.result[4, ] <- posterior.summary(sims.matrix.Kan[, "lambda"])
each.PM.lambda.result[5, ] <- posterior.summary(sims.matrix.Noda[, "lambda"])
each.PM.lambda.result[6, ] <- posterior.summary(sims.matrix.Abe2[, "lambda"])
round(each.PM.lambda.result, 3)

png("Figure_S4.png", width = 100, height = 60, units = "mm", pointsize = 8, res = 800)
par(mar = c(5, 2, 2, 2))
dotchart(rev(each.PM.lambda.result[, 1]), 
         labels = c("Abe (2)", "Noda", "Kan", "Fukuda", "Abe (1)", "Koizumi"), 
         pch = 19, xlim = c(-0.05, 0.1), xlab = "Effects (Point)", ylab = "")
points(each.PM.lambda.result[, 1], c(6:1), pch = 19, col = color[12])
segments(each.PM.lambda.result[, 2], c(6:1), each.PM.lambda.result[, 3], c(6:1), col = color[12])
abline(v = 0, lty = 3)
dev.off()

#### test for politically motivated sse of reshuffle cues (online Appendix H) ####
# preparing data for whether left- and right-leaning firms provided reshuffle cues
reshuffle.cue.pattern.binary <- matrix(NA, 11, 13)
reshuffle.cue.pattern.2 <- reshuffle.cue.pattern[1:11, 1:13]
reshuffle.cue.pattern.binary[reshuffle.cue.pattern.2 == "1" | 
                               reshuffle.cue.pattern.2 == "1*" | 
                               reshuffle.cue.pattern.2 == "1_" | 
                               reshuffle.cue.pattern.2 == "1_*"] <- 1
reshuffle.cue.pattern.binary[reshuffle.cue.pattern.2 == "0"] <- 0
partisan.cue.data <- data.frame(cue = c(reshuffle.cue.pattern.binary[3, ],  # Yomiuri
                                        reshuffle.cue.pattern.binary[6, ],  # Sankei-FNN
                                        reshuffle.cue.pattern.binary[7, ],  # Nikkei
                                        reshuffle.cue.pattern.binary[11, ],  # NNN
                                        reshuffle.cue.pattern.binary[4, ],  # Asahi
                                        reshuffle.cue.pattern.binary[5, ],  # Mainichi
                                        reshuffle.cue.pattern.binary[10, ]),  # ANN
                                right.leaning = c(rep(1, 13 * 4), rep(0, 13 * 3)), 
                                DPJ.regime = rep(c(rep(0, 6), rep(1, 5), rep(0, 2)), 7), 
                                after.2009 = rep(c(rep(0, 6), rep(1, 7)), 7), 
                                house = rep(c(1:7), each = 13), reshuffle = rep(c(1:13), 7))

## ideological position of polling firms and the use of reshuffle cues (Table S7)
# overall period (Panel A)
table(partisan.cue.data$DPJ.regime, partisan.cue.data$cue, partisan.cue.data$right.leaning)

# after 2009 (Panel B)
table(partisan.cue.data$DPJ.regime[partisan.cue.data$after.2009 == 1], 
      partisan.cue.data$cue[partisan.cue.data$after.2009 == 1], 
      partisan.cue.data$right.leaning[partisan.cue.data$after.2009 == 1])

## linear model for the use of reshuffle cues (Table S8)
# overall period (first column)
partisan.cue.result.1 <- lm(cue ~ right.leaning * DPJ.regime, data = partisan.cue.data)
partisan.cue.result.1.vcov <- cluster.vcov(partisan.cue.result.1, cluster = ~ house + reshuffle)
round(coeftest(partisan.cue.result.1, vcov. = partisan.cue.result.1.vcov), 3)
length(partisan.cue.result.1$residuals)

# after 2009 (second column)
partisan.cue.result.2 <- lm(cue ~ right.leaning * DPJ.regime, 
                            data = subset(partisan.cue.data, after.2009 == 1))
partisan.cue.result.2.vcov <- cluster.vcov(partisan.cue.result.2, cluster = ~ house + reshuffle)
round(coeftest(partisan.cue.result.2, vcov. = partisan.cue.result.2.vcov), 3)
length(partisan.cue.result.2$residuals)

#### Table S1 (estimated results of additional analyses) ####
## column (1)
print(coef.table(sims.matrix.inauguration[, 763:ncol(sims.matrix.inauguration)]), quote = FALSE)

## column (2)
print(coef.table(sims.matrix.self.selection[, 763:ncol(sims.matrix.self.selection)]), quote = FALSE)

## column (3)
print(coef.table(sims.matrix.varying[, 763:ncol(sims.matrix.varying)]), quote = FALSE)

## column (4)
print(coef.table(sims.matrix.two.types[, 763:ncol(sims.matrix.two.types)]), quote = FALSE)

#### Table S4 (AR(1) and AR(2) models) ####
## AR(1) model
print(coef.table(sims.matrix.AR1[, 763:ncol(sims.matrix.AR1)]), quote = FALSE)

## AR(2) model
print(coef.table(sims.matrix.AR2[, 763:ncol(sims.matrix.AR2)]), quote = FALSE)

#### Table S5 (estimated results of robustness checks) ####
## column (1)
print(coef.table(sims.matrix.wo.Kyodo[, 763:ncol(sims.matrix.wo.Kyodo)]), quote = FALSE)

## column (2)
print(coef.table(sims.matrix.party[, 763:ncol(sims.matrix.party)]), quote = FALSE)
# delta'
print(coef.table(sims.matrix.party[, 777:787] - sims.matrix.party[, 766:776]), quote = FALSE)

## column (3)
print(coef.table(sims.matrix.Nikkei1[, 763:ncol(sims.matrix.Nikkei1)]), quote = FALSE)

## column (4)
print(coef.table(sims.matrix.Nikkei2[, 763:ncol(sims.matrix.Nikkei2)]), quote = FALSE)

## column (5)
print(coef.table(sims.matrix.NNN[, 763:ncol(sims.matrix.NNN)]), quote = FALSE)

## column (6)
print(coef.table(sims.matrix.wo.kappa[, 763:ncol(sims.matrix.wo.kappa)]), quote = FALSE)

#### Table S6 (analyses for a signle PM) ####
## column (1)
print(coef.table(sims.matrix.Koizumi[, (n.weeks[1] + 1):ncol(sims.matrix.Koizumi)]), quote = FALSE)

## column (2)
print(coef.table(sims.matrix.Abe1[, (n.weeks[2] + 1):ncol(sims.matrix.Abe1)]), quote = FALSE)

## column (3)
print(coef.table(sims.matrix.Fukuda[, (n.weeks[3] + 1):ncol(sims.matrix.Fukuda)]), quote = FALSE)

## column (4)
print(coef.table(sims.matrix.Kan[, (n.weeks[6] + 1):ncol(sims.matrix.Kan)]), quote = FALSE)

## column (5)
print(coef.table(sims.matrix.Noda[, (n.weeks[7] + 1):ncol(sims.matrix.Noda)]), quote = FALSE)

## column (6)
print(coef.table(sims.matrix.Abe2[, (n.weeks[8] + 1):ncol(sims.matrix.Abe2)]), quote = FALSE)